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Abstract 

In  this  paper  we  analyse  a  sub-domain  residual  error  estimator  for  finite 
element  approximations  of  elliptic  problems.  It  is  obtained  by  solving  local 
problems  on  patches  of  elements  in  weighted  spaces  and  provides  for  an  upper 
bound  on  the  energy  norm  of  the  error  when  the  local  problems  are  solved  in 
sufficiently  enriched  discrete  spaces.  A  guaranteed  lower  bound  on  the  error 
is  also  derived  by  a  simple  postprocess  of  the  solutions  to  the  local  problems. 
Numerical  tests  show  very  good  effectivity  indices  for  both  the  upper  and  lower 
bounds  and  a  strong  reliability  of  this  estimator  even  for  coarse  meshes. 


1  Introduction 

The  use  of  residual-based  methods  for  a  posteriori  estimation  of  errors  in  finite  el¬ 
ement  approximations  of  elliptic  boundary-value  problems  has  become  important 
in  applications  as  methods  for  assessing  quality  of  numerical  solutions  and  provid¬ 
ing  a  basis  for  adaptive  meshing.  A  survey  of  such  methods  can  be  found  in  the 
monograph  of  Ainsworth  and  Oden  [1],  Among  early  residual  methods  was  the 
subdomain-residual  method  of  Babuska  and  Rheinboldt  [3],  which  involves  solving 
local  problems  on  patches  of  elements  for  which  the  local  residuals  appear  as  data. 

One  issue  that  frequently  complicates  residual-based  methods  is  the  treatment  of 
boundary  conditions  on  elements  or  patches  in  formulating  the  local  problems  for  er¬ 
ror  measures.  In  the  equilibrated  element  residual  methods  of  Ladeveze  and  Leguil- 
lon  [11]  and  Ainsworth  and  Oden  [2],  Neumann  conditions  are  used  and  the  residual 
fluxes  are  equilibrated  for  each  element.  This  process  can  lead  to  good  error  estima¬ 
tors,  but  is  generally  expensive  and  rarely  used  in  three-dimensional  applications. 
The  use  of  homogeneous  Dirichlet  boundary  conditions,  however,  lead  to  global  lower 
bounds,  and  often  results  in  estimators  that  do  not  deliver  acceptable  accuracy  for 
use  in  adaptive  mesh  strategies. 

More  recently,  variations  in  the  subdomain-residual  method  have  been  proposed  by 
Carstensen  and  Funken  [6]  and  Morin,  Nochetto  and  Siebert  [12]  (see  also  Datta 
[8]),  that  employ  a  partition  of  unity  generated  by  the  usual  finite  element  basis 
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functions  for  triangular  and  quadrilateral  meshes.  By  weighting  the  local  bilinear 
forms  using  such  partitions  of  unity,  local  problems  for  error  indicators  are  formed 
on  the  support  (patch-domains)  of  the  basis  functions,  and  the  weight  functions 
naturally  vanish  on  the  patch  boundaries.  This  provides  a  natural  and  convenient 
framework  for  overcoming  problems  of  assigning  boundary  averages  or  equilibrating 
residual  fluxes  on  the  boundaries. 

In  the  present  work,  we  present  a  detailed  analysis  of  a  variant  of  the  method  of 
Morin  et  al.  [12]  and  demonstrate  its  performance  on  several  model  problems  on 
one-  and  two-dimensional  domains.  We  develop  a  new  computable  error  estimator 
and  derive  both  upper  and  lower  bounds  on  the  global  error  measured  in  the  energy 
norm.  Our  numerical  experiments  suggest  the  method  is  easily  implemented  for 
elliptic  problems  and  can  yield  excellent  bounds  on  the  global  error. 


2  Model  problem  and  preliminaries 

Let  O  be  an  open  bounded  polygonal  domain  in  Md,  d  =  1  or  2,  with  boundary  30 
and  let  u:  O  — *  R  be  the  solution  of  the  problem 


—Art  +  cu  =  f 
u  =  0 

where  /  £  L2(0)  and  c  is  a  positive  constant, 
reads: 

Find  u  £  Hq  (O)  such  that  /  (Vu  ■  Vv  +  cuv )  =  /  fv  V  v  £  Hq  (O)  (2) 

J  Q  J  Q 

where  ILq(O)  is  the  usual  Sobolev  space  of  functions  which  are  in  L2(0)  with  square 
integrable  derivatives  and  that  vanish  on  30.  The  bilinear  form  -£>(•,•)  defined  on 
^(O)  x  H'( O)  by 

B(u,v)  =  /  (Vu  ■  Vv  +  cuv)  (3) 

Jn 

is  symmetric,  continuous  and  coercive.  It  thus  defines  an  inner  product  in  H d(0) 
and  induces  the  so-called  energy  norm  |||u|||  =  \jB(v,v).  Introducing  the  bounded 
linear  form  F(-)  defined  on  ILg(O)  as 

F(v)  =  [  fv, 

Jn 

the  Lax-Milgram  Theorem  (see  e.g.  [10])  allows  us  to  establish  that  the  problem: 


in  O 
on  30 


(1) 


The  weak  formulation  of  Problem  (1) 


Find  u  £  iLg(O)  such  that  B(u,v)  =  F(v),  V  v  £  Hq(Q) 


(4) 
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is  well-posed  in  the  sense  that  it  admits  a  unique  solution  u  G  (SI),  which  de¬ 
pends  continuously  on  the  data.  Moreover,  u  is  the  solution  of  (1)  in  the  sense  of 
distributions. 

Let  Th  be  a  regular  finite  element  partition  of  made  up  of  triangular  or  quadrilat¬ 
eral  elements  (in  2D).  We  denote  by  h k  the  diameter  of  each  element  K  of  the  mesh 
Th  and  by  pk  the  diameter  of  the  largest  circle  inscribed  in  K .  if  K  is  a  triangle. 
A  similar  definition  holds  if  K  is  a  quadrilateral  [1,  7].  The  mesh  is  assumed  to  be 
regular,  i.e.  there  exists  a  positive  constant  k  >  0  such  that: 

Kk  =  — —  <  K  V  K  G  Th-  (5) 

Pk 

We  then  introduce  the  finite  dimensional  space  C  Hq(FI)  that  is  the  usual  space 
of  Pfc  (for  meshes  of  triangles)  or  (for  meshes  of  quadrilaterals)  continuous  finite 
elements  (see  e.g.  [7]  or  [10]  for  a  precise  definition).  The  finite  element  approxima¬ 
tion  of  problem  (4)  reads: 

Find  Uh  G  such  that  B(uh,v)  =  F(v)  V  v  G  Vj£.  (6) 

Denoting  the  approximation  error  in  Uh  by  e  =  u—Uh,  we  have  the  following  equation 
for  the  error 

B(e,v)  =  B(u  —  Uh,v)  =  F(v)  —  B(uh,v)  =  TZ(v)  V»G  Hq(Q).  (7) 

where  77(.)  is  called  the  residual  functional.  Note  that  the  well-known  orthogonality 
property  reads: 

K(v)  =  B{e,  v)  =  0  VveV%.  (8) 


Lemma  1  With  above  assumptions  and  definitions,  the  following  equality  holds: 

\K(v)\  _ 


|e|||  =  sup  =  \\K\\i 

veHi(n)  lll^lll 


(9) 


where  ||77||*  is  the  norm  of  the  residual  in  Ft  1(D). 

Proof:  Using  (7),  we  have 

m  m  B(e,e)  77(e)  ^  |77(u)| 

lllelll  =  =  Til — TiT  ^  SUP  TiTTiT- 

lllelll  lllelll  veHi(Q)  lllelll 

Furthermore,  from  the  continuity  of  £?(•,•),  |77(u)|  =  \B{e,v)\  <  |||e|||  |||u|||  for  all 
v  G  Ft o(D),  which  establishes  the  assertion.  □ 


Note  that  the  equality  in  relation  (9)  holds  because  the  bilinear  form  B(-,-)  is 
positive  definite  and  symmetric.  In  more  general  cases,  we  should  expect  to  have 
only  an  equivalence  between  the  norm  of  the  residual  and  the  norm  of  the  error. 
Relation  (9)  will  be  useful  in  deriving  a  lower  bound  on  the  error. 
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Figure  1:  Representation  of  the  function  fa  on  ojt  for  triangular  (left)  and  quadran¬ 
gular  (right)  meshes. 

3  Upp  er  bound  estimate  on  the  error 

Let  fa  be  the  set  of  Lagrangian  piecewise  linear  (for  triangles)  or  piecewise  bilinear 
(for  quadrilaterals)  basis  functions.  The  support  of  each  <j>i  is  denoted  by  coi  and 
will  be  referred  to  as  the  patch  of  elements  connected  to  node  i  of  the  mesh  (see 
Fig.  1  for  a  representation  of  fa  and  u in  two  dimensions).  Subsequently,  an  interior 
patch  uji  will  define  a  patch  for  which  node  i  is  not  located  on  the  boundary  <9P. 
Similarly,  a  boundary  patch  Ui  is  a  patch  for  which  node  i  is  on  dCl. 

We  will  use  the  fact  that  the  functions  fa  form  a  so-called  partition  of  unity 

^  fa(x)  =  1  VrGll  (10) 

l<i<N 

where  N  denotes  the  total  number  of  nodes  of  r^.  Furthermore,  let  hi  =  maxjf£w;  hx 
and  pi  =  min^e^j  Pk ■  The  regularity  of  the  partition  t/j  is  inherited  by  the  patches 
ujt  and  we  have  (see  Theorem  6.1  in  [1]): 

3  Co  >  0:  Ki  =  —  <  Co  k  V  i  =  1  (11) 

Pi 


Note  that  upon  inserting  (10)  into  (7),  we  have 

B(e,v)=K(v  Y,  (t>i)=  J2  VvG^(fi).  (12) 

l<i<AT  l<i<N 

Let  us  introduce,  now,  on  each  patch  c Uj,  the  bilinear  form 

v)=  fa{Vu  ■  Vu  +  cuv )  (13) 

J  UJi 

with  associated  norm  HHH^  =  y/B^v,  v),  and  define  the  space  for  each  i  = 
1 ,N  as 

4*1 


=  W{LOi) 
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where  for  an  interior  patch  uf 


yv(wi)  =  c1^) 


and  for  a  boundary  patch  ujf. 

W(uJi)  ={t)£  v  =  0  on  du>i  n  90}. 

Next,  we  will  formulate  and  solve  the  following  N  local  problems,  over  each  patch 
u>i,  for  local  error  functions  Q  G  W(uf) 


B4>i{QiA)  =  1  <  i  <  N. 

We  define  the  global  error  estimator 


(14) 


e  = 


1  E 

l<i<JV 


(15) 


We  show  in  the  following  that  i  provides  for  an  upper  bound  on  the  error  |||e 


Remark  1  The  left-hand  side  and  right-hand  side  of  (14)  can  be  explicitly  written 
as: 


B<t>MiA)  =  /  (0jVO  •  +  ccfiCiip), 

TZ(^<f>i)=  f  ffxfi-  f  (Vuft-V(#j)  +  cMft#i). 
J  Ui  J  LOi 


(16) 


Whenever  uoi  is  an  interior  patch,  no  boundary  conditions  are  prescribed  on  duii  in 
(14)-  Thus,  the  local  problems  (14)  are  of  11 Neumann ”  type.  In  the  case  c  =  0,  the 
solution  is  defined  up  to  a  constant  and  we  should  expect  to  have  a  compatibility 
condition  on  the  right-hand  side.  Indeed,  thanks  to  the  Galerkin  orthogonality  (8), 
it  happens  that  for  if  a  constant  function  on  LOi 


KWh)  =  if>K(<f>i)  =  0.  (17) 

Thus,  the  compatibility  condition  is  always  satisfied.  In  order  to  fix  the  constant,  we 
seek  a  solution  of  (14)  in  the  space 

W(ui)  =  {v  eW(ui),  [  v<t>i  =  0}  (18) 

J  LOi 


instead  ofW(ui). 

Remark  2  As  written  in  (16),  the  residual  does  not  contain  any  integrals  of  the 
fluxes  along  the  element  interfaces.  As  presented,  it  is  shown  that  it  is  not  necessary 
to  perform  the  integration  by  parts  as  most  authors  do  [6,  12]. 
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In  Lemma  3,  we  will  prove  the  well-posedness  of  the  local  problems  (14).  Toward  this 
goal,  we  state  Lemma  2  below.  Its  proof  has  been  given  in  [12]  for  a  triangular  mesh 
and  is  provided  in  the  Appendix  in  the  case  of  quadrangular  meshes.  We  remark  that 
the  proof  varies  substantially  when  dealing  with  triangles  or  quadrilaterals.  Indeed, 
the  weighting  function  fa  is  linear  for  the  first  type  of  elements.  For  quadrilateral 
elements,  fa  is  bilinear  which  means  that  the  function  can  be  quadratic  along  certain 
directions  in  the  elements,  making  the  proof  a  little  more  complex. 

Lemma  2  (Weighted  Poincare  inequality)  Let  hi  be  the  maximum  size  of  the 
elements  in  the  patch  Wj.  For  any  function  fa  £  W{uf)  that  satisfies  the  condition 
f  fafa  =  0  whenever  uii  is  an  interior  patch,  there  exists  a  constant  C,  independent 
of  fa,  but  dependent  of  n,  such  that 

1 1  Ci  1 1 L2  (u)i )  <  ChiWVfaW^,  (19) 


where  ||  •  ||0i  =  ||  •  fa1/2 ||£2(wi). 

Lemma  3  Problem  (14)  is  well  posed  for  any  given  hi  >  0. 


Proof:  We  first  proof  this  lemma  for  an  interior  patch  tu*.  The  bilinear  form  ( • ,  •) 
is  an  inner  product  on  W(u>i)  and  thus  is  continuous  and  coercive.  The  assertion 
then  follows  from  the  Riesz  representation  theorem  if  we  can  show  that  the  functional 
lZ(ipfa)  is  bounded  in  Wfaif)  with  respect  to  the  norm  |||V;|||</>i- 

O 

Let  us  first  remark  that  if  c  =  0,  all  functions  fa  £  Wfauf)  satisfy  the  condition 
f  fafa  =  0.  On  the  other  hand,  if  c  fa  0,  thanks  to  the  Galerkin  orthogonality  (8), 
we  have  IZ(fa)  =  0  and,  for  all  fa  £  Wfajf), 

f  fafa  \ 

=W«). 

u  &  J 


IZfaffa)  =  u 


where  fa  satisfies  now  the  condition  fafai  =  0.  Hence, 


\R(fafa)\  =  \K(fafa)\  = 

ffafa- 


[  ffafa  ~  [ 

J  LUo  j  (jJ 


(V%  •  V(V^)  +  cuh^(j>i) 


< 


faJuh  ■  Vfa  +  cuhfafa 


J 

J  U>i 


faVuh  ■  Vfa 


<  WfbMUi  +  W^uhUiWVfaWfr  +  c||«fc||^ 


/ 

J  UJn 


faVuh  ■  Vfa 


(20) 


Let  us  remark  now  that 


(21) 
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Indeed,  if  we  set  g  =  Ju_  ifrfa/ fu.  fa,  we  have 


llffl 


2 

4*1 


(Li Mi)2  K  MlliiLjfo) 

Li  &  Li  & 


IIV'll 


2 


and  clearly, 

Wi’hi  <  M k  +  hhi  <  2I 

Thus,  noticing  that  X7ip  =  X7ip  and  using  the  previous  inequality,  we  have 


mfa)\  <  2||/||^||^||^  +  ||Vuh|k||V^|k  +2cK|k 


/ 

J  UJi 


tpVuh  ■  V fa 


We  now  need  to  bound  the  term 


f  fas/uh  ■  Vfa  .  We  have  by  Holder’s  inequality: 


'ipVuh  ■  Vfa 


lUJi 


<  ||V0i|Uoc(wiI  f|V«h||L2(wi)  IMI La(wi) 


Moreover,  a  well-known  property  of  Lagrangian  basis  functions  allows  us  to  state 
that 

\\Vfa\\L°°(in  <  C i —  MK  G  Ui  for  K  quadrilaterals, 

Pi< 

C2 

\\V fall l°° ( k)  <  —  Vit  G  Wj  for  K  triangles. 

PK 


In  both  cases, 


C' 

||V0i|| L°°(ui)  =  max  ||V&||L=o(tf)  <  — 

KGcJi  Pi 


(22) 


where  p,,  =  min^e(Ji(/9x)  and  C'  might  depend  on  the  regularity  parameter  k  of  the 
partition  t/(i. 


Thus, 


/ 

LJ,- 


ijj'VuhVfa 


C' 

-  “llVuklU2(wi) 

Pi 

C'Ch 

<  — —  UVufcll^^)  IIV^IU 

Pi 

C'Ch 

<  — -illVufcll^^)  IIVV-IU 

Pi 


where  we  have  used  the  weighted  Poincare  inequality  (19). 
Then,  we  have 


IW*)|<2||/|UJ|V’IUi  +  ||V^||^J|VV-|Ui+2c||^||0J|^||^ 

C'Ch  ■ 

+  — — ^iiVtifciiL^)  IIV^IU 
Pi 
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and  noticing  that  ||V^||^  <  HIVIIU  and  ||V’||</>i  <  c  ^WV’IIU,  we  finally  obtain 
/  o  d'Ch  \ 

I^WOI  <  y-^Wfhi  +  \\^uhUi  +  2Vc  IKIU  +  —  - ||v^t||L2K)j  Hi^llk 

(23) 

The  proof  is  thus  complete  since  all  the  quantities  appearing  on  the  right-hand  side 
are  bounded.  The  same  proof  is  also  valid  for  a  boundary  patch,  without  having  to 
introduce  the  function  iJj.  □ 

We  now  derive  an  exact  upper  bound  for  the  error.  Using  the  previous  definitions 
and  the  Cauchy-Schwarz  inequality,  we  have  for  all  v  £  (O) 


B{e,v)  =  Y 

1  <i<N 


=  Y 

1  <i<N 


1  <i<N 


< 


Taking  v  =  e  in  the  previous  inequality,  noticing  that  Xa <i<NB4>i(viv)  =  B(v,v) 
for  any  v  in  Hq(£1),  and  using  the  definition  of  i  in  (15),  we  obtain  the  following 
guaranteed  upper  bound  on  the  error,  i.e. 


=  \J B(e,  e)  <  e. 


However,  the  quantity  e  cannot  be  computed  exactly  because  the  local  problems 
(14)  defined  on  each  patch  for  Q  are  of  infinite  dimension.  We  show  in  the  next 
section  how  to  approximate  these  problems  in  order  to  derive  an  a  posteriori  error 
estimator. 


4  The  a  posteriori  error  estimator 

Let  Pk+p(uJi)  denote  the  finite  element  space  of  piecewise  polynomials  of  degree  k+p 
on  the  patch  ul  and  let  Vk+p(iOi)  =  Pk+p{uJi)  fl  W{u)j).  Here,  p  represents  the  extra 
degree  with  respect  to  the  degree  k  used  to  compute  the  finite  element  solution  u /, . 

On  each  patch  coi,  we  now  solve  for  r/7;  £  Vk+p{iOi)  as  the  solution  to  the  problem 

=  ^(#»)  V  V’  G  Vk+p(cOi)  (25) 
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and  define  the  global  error  estimator 

e=  /  Y  £2i  =  J  H  (26) 

y  1  <i<N  y  l<i<N 

The  quantities  £*  =  i ^/n^Jrjifrji)  represent  the  local  contribution  from  each  patch 
uit  to  the  global  error  estimator  e. 


Remark  3  Although  the  residual  vanishes  on  we  can  actually  choose  p  =  0  in 
(25)  because  the  term  IZ(ip(pi)  does  not  necessarily  vanish  for  if  £  Vk{iOi)  due  to 
the  presence  of  the  weighting  function  (pi.  However,  we  expect  to  get  more  accurate 
estimates  as  p  is  increased.  Indeed  we  can  verify  from  Problem  (25)  that  the  error 
estimator  e  should  increase  as  the  value  of  p  increases  since: 


= 


sup 

veVk+p(u>i) 


\T^{v4>i)\ 
IIMIk  ' 


(27) 


Lemma  4  Let  Yl)l+p(Ll)  be  the  global  finite  element  space  of  piecewise  polynomials 
of  degree  k  +  p  defined  on  and  vanishing  on  dPl.  Then  we  have 

|||e|||<e  +  2  inf  |||it  —  u|||.  (28) 

v£Vlh+p(0.) 


Proof:  On  one  hand,  we  have  for  all  v  £  vj)+p(fi), 

Ik  -  uh\\\2  =  B{y  -uh,v-  uh ) 

=  B(e,v  -  uh)  +  B{v  -  u,v  -  uh)  (29) 

<  B(e,v  —  Uh)  +  \  \\v  —  u\\\  |||u  —  Uh\\\- 

On  the  other  hand,  using  (7)  and  the  Cauchy-Schwarz  inequality,  we  have 

B(e,  v  -  uh)  =  -  uh) 

=  V  K{(v  -  uh)4>i) 

1  <i<N  (30) 

=  Y  n({v  -  UlO^Cpi) 

l<i<N 

where  (v—Uh)\u,i  denotes  the  restriction  of  ( v—Uh )  on  Wj.  Since  (v—Uh) \UJi  £  Vk+p(iOi) 
and  using  (25),  we  may  further  write 

B(e,v-uh)=  Y  -Uh)^.) 

l<i<N 


< 


l<i<N 


BfaiVuVi)  Y  Bfciv  —  uh,  v  —  uh) 


(31) 


l<i<N 
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With  (29)  and  (31),  we  obtain  |||u  —  <  e  +  |||u  —  n|||  and  |||e|||  <  |||it  —  v|||  + 

| ||u  —  Uh\ ||  <  £  +  2|||w  —  u|||.  Since  the  last  inequality  holds  for  all  v  £  vj)+p(0), 
inequality  (28)  is  thus  established.  □ 

Inequality  (28)  shows  that  the  estimator  e  does  not  necessarily  provide  for  an  upper 

bound  on  the  error.  However,  the  quantity  2  inf  \\\u  —  u|||  should  decrease  as  p 

veVk+p(Q) 

increases,  so  that  e  should  be  an  upper  bound  for  high  enough  p. 

As  a  final  result,  we  show  that  the  estimator  e  is  equivalent  to  the  energy  norm  of 
the  error. 

Lemma  5  There  exists  a  constant  C  >  1,  that  depends  only  on  the  regularity  of  the 
mesh  and  on  the  constant  C  appearing  in  Lemma  2,  but  otherwise  independent  of  h 
and  p,  such  that 

e  <  C'|||e|||.  (32) 


Proof:  For  any  i  =  1, . . . ,  N,  the  function  r)i<pi  is  a  piecewise  polynomial  function 
that  vanishes  on  the  boundary  diOi.  Let  us  denote  by  rn(f>l  the  extension  of  rg/fi  by 
zero  over  H  \ uj{.  Clearly,  rgfi  belongs  to  Hq(Q)  and  can  be  taken  as  a  test  function 
in  the  continuous  problem  (4).  Moreover,  its  gradient  will  also  vanish  outside  tUj. 
We  thus  have, 


Hence, 


fVi<t>i  =  /  fVi^i  =  B(u,  rjicfi)  =  /  (Vu  •  V(rn(j)i)  +  cwpcfi) 

WhiWlli  =  i<t>i)  =  /  fVi&i  -  /  (Vttfc  •  V(rfi(f>i)  +  cuhr}i4>i) 
J  LOi  J  UJi 

=  /  [V (u  -  uh)  ■  V (rjiffi)  +  c(u  -  uh)rii(j)i\ 

J  LOi 

=  B<j>i(u-uh,rn)+  V{u  -  uh)  -Vfarii. 

J  LOi 


(33) 


Let  us  remark  now  that  the  local  solution  r/,;  belongs  to  the  space  W(ui)  also  in  the 
case  c  >  0.  Indeed,  by  taking  a  constant  test  function  if  =  c  in  (25)  we  have 

c)  =  c  /  rjicfa  =  7 l(c<j>i)  =  c  7 Z{<j>i)  =  0, 

J  UJi 

owing  to  the  Galerkin  orthogonality.  The  local  functions  r/i  then  satisfy  the  weighted 
Poincare  inequality, 

H?7i|lL2(^)  <  ChiWVrjiWfc. 

We  can  proceed  in  a  similar  way  as  in  the  proof  of  Lemma  3  to  bound  the  last  term 
in  (33),  thus  obtaining 


«/*IIUIII^IIk  +  C'||V(«-«h)||L2(h,l)[|Vr7i||0i, 
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where  the  constant  C  depends  on  the  constant  appearing  in  the  weighted  Poincare 
inequality  and  on  the  regularity  parameter  k  of  the  mesh.  Summing  over  all  the 
patches  u>i  and  using  the  discrete  Cauchy-Schwarz  inequality,  we  have 


e2<  J  W\u~uh\ 

That  is: 


l<i<N 


J2  £t+cJ  J2  iiv(«-«h)iii2(w0  /  £l 


l<i<N 


l<i<N 


l<i<N 


£  <  \\\u  -  uh\\\  +  C  fiV(ii-uh)|| 


2 

L2(ui)' 


(34) 


1<  i<N 


For  each  element  K  of  the  mesh  77,, ,  let  us  denote  Nk  =  {i  ■  A'  C  wj,  A k  =  dim  Nk 
and  A  =  rna x^er/i  A k-  Then,  we  have 


J2  WV(U~Uh)\\h(oJi)=  \\V(u-uh)\\2L2{K) 

1  <i<N  KCuJi 

=  XkWV(U~  Uh)\\h(K)  <  A|||«-«ft| 


l<j<JV 


(35) 


KCth 


and  the  final  result  (32)  is  achieved  with  (7  =  1  +  Cy/\. 


□ 


5  Recovery  of  a  global  lower  bound  on  the  error 


Relation  (9)  shows  that  for  any  function  z  €  Hq(Q),  the  quantity  |7£(z)|/||M|| 
provides  a  lower  bound  of  the  error.  Here  we  choose  to  construct  the  function  z  as 

z=  Vifa, 

l<i<N 


r]i  being  the  solution  of  problems  (25).  The  function  z  is  continuous  since  the  basis 
functions  fa  vanish  at  the  boundaries  of  all  patches;  hence,  z  is  in  H g(H).  This 
choice  is  motivated  by  the  fact  that  the  functions  77*  are  already  available  upon 
solving  the  local  problems  (25). 

Furthermore,  we  note  from  (25)  that: 

K(z)  =  K(  ^2  rli(t>i)=  B(rli,Vi)  =  £2- 

l<i<N  l<i<N  l<i<N 

The  lower  bound  on  the  error  then  reduces  to: 

T]iow  —  m^iii  •  (36) 

In  other  words,  since  e2  is  already  known,  it  suffices  to  compute  the  norm  |||^|||  = 
II  Xyi<i<iv  Vifiil  1 1  t°  derive  a  lower  bound  on  the  error. 
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Remark  4  As  observed  in  Remark  1,  in  the  case  c  =  0,  the  solutions  rji  of  prob¬ 
lem  (25)  on  an  interior  patch  uii  are  obviously  defined  up  to  an  arbitrary  constant. 
We  note  that  the  upper  bound 


[  iVr/*i2^ 

1  <i<N 

does  not  depend  on  this  “ arbitrary ”  constant.  On  the  other  hand,  for  the  lower 
bound  rjiow,  the  quantity 


does  depend  on  the  constant.  Following  the  suggestion  in  Remark  1,  we  choose  the 
constant  Ci  as: 


Ci  = 


(37) 


O 

In  this  manner,  rg  —  Ci  G  Vk+P(u jfi)  n  W(ui)  and  the  function  z  is  defined  as: 


Z  =  ini  ^  Ci)(fi. 

1  <i<N 


(38) 


6  A  strategy  for  mesh  adaptation 

In  addition  to  accuracy  assessment,  a  posteriori  error  estimators  can  be  used  as 
refinement  indicators  for  mesh  adaptation  in  order  to  reduce  the  numerical  error. 
A  simple  strategy  consists  of  defining  local  contributions  to  the  global  error  and 
refining  the  elements  which  exhibit  large  source  of  error. 

For  the  subdomain-based  error  estimator  proposed  here,  there  are  actually  two  ways 
for  splitting  the  global  error  into  sums  of  local  contributions:  namely,  we  can  define 
patch- wise  or  element-wise  contributions.  From  (26),  we  have 


where  e*  =  y/B(f>i(r]i,r]i)  represents  the  contribution  of  the  patch  ut  to  the  error 
estimate  e.  Similarly,  we  can  decompose  e  into  a  sum  of  element-wise  contributions. 
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Indeed, 


£  = 


J  Y  £i  =  \PY  [  (\VVi\2 +  cv?)<i>i 

y  i<i<iv  V  i<i<NJuJi 


N 


Y  (  Y 


KCrh  \i£NK 


/£■ 

K<Zrh 


-K 


where,  now,  the  quantity  £k  represents  the  element-wise  contribution.  The  pro¬ 
cedure  can  naturally  be  applied  to  the  lower  bound  estimator  as  well;  in  fact,  we 
have: 


Vlow  — 


E 


'-K 


\<K<Ne'''Z' 


Y  £k ■ 


low 


l<K<Ne 


with  £k,Iow  =  £a'/III2:III  local  contribution  per  element. 


In  the  numerical  experiments  shown  in  the  next  sections,  we  only  consider  refine¬ 
ments  based  on  the  element-wise  contributions  (in  order  to  be  able  to  compare  our 
mesh  adaptation  procedure  with  the  element  residual  method).  Namely,  an  element 
K  C  Th  is  refined  if  either 


£a 

max  (eft-) 
Kcrh 


—  Cadap 


or 


£ K,low 

max  (eK  low) 

kCTh 


E  Cadap 


(39) 


where  Cadap  is  a  user-prescribed  constant  parameter  ranging  from  zero  to  one.  In 
the  numerical  tests  presented  in  section  8,  we  have  chosen  Cadap  =  0.5. 


7  Numerical  experiments  in  ID 

The  primary  objective  in  these  examples  is  to  assess  the  performance  of  the  upper 
and  lower  bound  estimates  on  the  numerical  error.  In  all  experiments,  the  quality  of 
these  estimates  will  be  measured  in  terms  of  the  effectivity  index,  namely  the  ratio 
between  the  error  estimate  and  the  exact  error  in  the  energy  norm  1 1 1  •  1 1 1 . 

We  first  consider  the  two-point  boundary-value  problem: 

— u"  +  cu  =  /  in  (0, 1) 

u  =  0  at  x  =  0  (40) 

u  =  0  at  x  =  1 

where  the  constant  c  will  be  chosen  as  either  one  or  zero. 
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7.1  Problem  with  smooth  solution 

In  our  first  example,  we  select  the  load  /  such  that  the  exact  solution  (see  Fig.  2) 
is  analytically  given  by: 

u(x)  =  ^(1  —  x)(l  —  3x)  —  2x(x  —  ^)(1  —  x)2^j 

The  finite  element  approximations  are  computed  using  polynomial  degrees  k  =1,  2 
or  3  and,  for  each  case,  we  calculate  the  upper  bound  e  and  lower  bound  piow  setting 
the  extra  degree  as  p  =  0,  1,  2  or  3.  The  effectivity  indices  are  shown  in  Figs.  3,  4 
and  5  respectively  on  sequences  of  uniformly  refined  meshes. 


exact  solution 


Figure  2:  Representation  of  the  exact  solution  u(x). 

We  remark  first  that  the  effectivity  indices  are  nearly  identical  for  p  =  1,  2  or  3.  In 
other  words,  this  shows  that  the  local  problems  rapidly  “saturate”  as  p  is  increased. 
More  surprisingly,  we  also  observe  that  most  of  the  results  are  acceptable  when 
taking  p  =  0.  This  is  an  important  observation  as  in  many  existing  finite  element 
codes,  there  is  no  access  to  polynomial  shape  functions  of  higher  degree  than  the 
ones  used  for  the  finite  element  solution  Uh •  This  aspect  of  the  error  estimator  may 
be  an  advantage  over  other  element  residual  type  methods  for  which  it  is  necessary 
to  take  p  >  0.  However,  if  we  choose  p  =  0,  we  are  not  guaranteed  that  the 
quantity  e  provides  for  an  upper  bound  on  the  error  (even  as  h  tends  to  zero).  In 
this  case,  more  sophisticated  techniques  may  be  envisaged,  that  aim  at  recovering 
a  guaranteed  upper  bound  by  estimating  the  error  ||Ci  —  ViWfa  in  the  solution  of  the 
local  problems,  as  well,  as  proposed  in  [4,  13].  We  will  not  investigate  this  approach 
in  the  present  work. 

Another  observation  deals  with  the  fact  that  the  error  bounds  barely  depend  on  the 
mesh  size  h  (at  least  when  the  mesh  is  not  coarse).  Furthermore,  the  effectivity 
index  of  the  upper  bound  decreases  as  k  goes  from  one  to  three.  In  fact,  we  find  an 
effectivity  index  close  to  1.22  for  k  =  1,  1.07  for  k  =  2  and  1.04  for  k  =  3.  On  the 
other  hand,  the  lower  bound  piow  always  provides  for  effectivity  indices  very  close 
to  one  (even  for  p  =  0)  even  for  coarse  meshes.. 
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Figure  3:  Global  effectivity  index  for  the  upper  (left)  and  lower  (right)  bounds  for 
the  case  k  =  1  and  c  =  0  (top),  c  =  1  (bottom). 


7.2  A  variant  of  the  error  estimator 

In  [12],  the  authors  proposed  to  solve  the  local  problems  (25)  using  homogeneous 
Dirichlet  boundary  conditions.  In  other  words,  defining  Vq  +p(cuj)  as  the  subspace 
of  Vk+P(u Ji)  such  as: 

V^+p(LOi)  ={»£  Vk+p(uJi),v  =  0  on  dui}, 

they  solve  for  rj j  6  Vq+p(l jj)  such  that 

Vi,v)  =  n(v&)  VuG  Vk0+p(uJi).  (41) 

Here  we  compare  the  error  estimators  (upper  and  lower  bounds)  when  using  either 
Dirichlet  or  Neumann  boundary  conditions.  For  the  case  c  =  1,  k  =  1,  and  taking 
p  =  0  or  1,  we  obtain  the  results  shown  in  Fig.  6. 

We  observe  that  the  effectivity  indices  for  both  the  upper  and  lower  bounds  converge 
at  a  slower  rate  with  respect  to  the  number  of  degrees  of  freedom  when  using  ho¬ 
mogeneous  Dirichlet  boundary  conditions.  We  conclude  that  the  subdomain-based 
error  estimation  method  is  more  reliable  when  we  do  not  impose  any  boundary 
conditions  on  the  patch  boundaries. 
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c  =  0,  k  =  2 
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Figure  4:  Global  effectivity  index  for  the  upper  (left)  and  lower  (right)  bounds  for 
the  case  k  =  2  and  c  =  0  (top),  c  =  1  (bottom). 


7.3  Comparison  with  other  error  estimators 

In  this  section,  we  compare  the  subdomain-based  error  estimator  (SRM  for  Subdo¬ 
main  Residual  Method)  with  other  methods:  namely,  the  ZZ  recovery-type  method 
[14,  15],  and  the  Element  Residual  Method  (ERM)  [5,  9].  We  choose  here  to 
solve  (40)  with  c  =  0.  Note  that  the  results  are  in  this  case  the  same  whether 
we  use  either  ERM  or  the  Element  Residual  Method  with  the  equilibration  tech¬ 
nique  (EqRM)  [11,  1], 

We  want  to  test  these  three  estimators  using  an  approximation  of  the  Green  function 
as  a  solution  of  (40).  We  recall  that  the  Green  function  is  the  solution  of  (40) 
(with  c  =  0)  when  the  load  /  is  the  Dirac  distribution.  The  Green  function  is 
particularly  interesting  because  it  possesses  a  discontinuity  in  the  first  derivative  at 
the  point  where  the  Dirac  distribution  is  located.  To  compute  a  solution  of  (40) 
which  approximates  the  Green  function,  we  actually  replace  the  Dirac  distribution 
located  at  x  =  xq  by  the  function  f  =  ke  defined  by: 
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Figure  5:  Global  effectivity  index  for  the  upper  (left)  and  lower  (right)  bounds  for 
the  case  k  =  3  and  c  =  0  (top),  c  =  1  (bottom). 
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Figure  6:  Global  effectivity  index  for  the  upper  and  lower  bounds,  Dirichlet  or 
Neumann  boundary  conditions. 


where  the  constant  C  is  chosen  to  satisfy 


/  kt{x)dx  =  1. 

Jo. 


(43) 


In  the  following,  we  take  e  =  0.1  for  which  the  corresponding  load  /  and  “exact 
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solution  u  are  plotted  in  Fig.  7.  In  this  case,  the  exact  solution  is  not  known  and 
is  approximated  by  an  overkilled  finite  element  solution  Uh  using  1500  degrees  of 
freedom. 


Dirac  approximation  with  epsilon  =  0.1  exact  solution 


Figure  7:  Representation  of  the  load  (left)  and  the  exact  solution  (right). 

In  order  to  compute  effectivity  indices,  we  also  calculate  an  accurate  value  of  the 
energy  norm  of  u  such  as: 


IIMH  ~  VWiW+WW 

where  e  is  the  error  in  u h  estimated  here  by  ERM.  The  “exact”  error  e  in  the 
finite  element  solution  Uh  will  then  be  calculated  by  the  formula  (obtained  from  the 
orthogonality  property  (8)): 

1 1  Mil2  =  IIMII2  -  IIMII2  ~  IIMII2  +  IIMH2  -  HMII2  (44) 


In  Fig.  (8),  we  present  the  different  results  obtained  with  the  three  error  estimators. 
In  each  case,  the  exact  error  is  estimated  using  (44).  We  take  k  =  1  for  the  degree 
of  the  finite  element  solution,  p  =  1  for  the  extra-degree  used  in  SRM  and  ERM. 
The  first  three  graphs  represent  the  distribution  of  the  error  estimate  on  the  domain 
[0, 1]  for  SRM,  ERM,  and  ZZ  using  30  elements.  We  observe  that  the  approximation 
error  is  concentrated  in  the  middle  elements,  where  the  main  variation  in  the  first 
derivative  of  the  exact  solution  actually  occurs.  The  last  graph  represents  the  be¬ 
havior  of  the  effectivity  index  for  these  three  estimators  when  the  number  of  degrees 
of  freedom  is  increased. 

Our  main  conclusion  is  that  the  ZZ  method  is  poorly  accurate  when  we  do  not  use 
enough  degrees  of  freedom  (less  than  50  degrees  of  freedom  in  our  example).  This 
is  due  to  the  fact  that  the  ZZ  method  is  based  on  a  recovery  of  the  gradient  of 
the  solution  Uh.  However,  the  two  other  estimators  seem  to  be  quite  stable.  The 
effectivity  index  with  SRM  stays  close  to  1.22,  which  is  the  value  we  obtained  in 
our  previous  experiments,  whereas  the  effectivity  index  using  ERM  remains  close  to 
unity. 
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Figure  8:  Distribution  of  the  error  estimate  with  ERM  (top  left),  SRM  (top  right) 
and  the  ZZ  (bottom  left).  Effectivity  indices  for  each  method  (bottom  right). 


8  Numerical  experiments  in  2D 

8.1  Performance  of  the  estimator 

We  now  investigate  the  performance  of  the  estimator  in  the  case  of  2D  problems. 
We  consider  the  squared  domain  Q  =  {(x,  y)  G  M2,  0  <  x  <  1,  0  <  y  <  1}  and  solve 
the  problem 

—Art  +  cu  =  f  in  il 

n  on  (45) 

u  =  (J  on  oil 

where  c  is  a  constant  chosen  as  either  0  or  1.  The  load  /  is  such  that  the  exact 
solution,  plotted  in  Fig.  9,  is  given  by: 

u(x )  =  0.0005  x2(l  —  x)2y2(l  —  y)2 e10^2  +y\  (46) 

The  following  figures  present  the  results  for  different  values  of  the  parameters:  the 
degree  of  the  finite  element  solution  is  chosen  as  k  =  1  or  2  whereas  the  extra  degree 
used  to  compute  the  upper  and  lower  bounds  can  take  the  value  p  =  0,  1,  2  or  3. 
As  in  the  ID  experiments,  we  show  the  effectivity  indices  of  the  error  estimates. 
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Figure  9:  Representation  of  the  exact  solution  u(x). 


The  results  are  shown  in  Figs.  10  and  11  for  k  =  1  and  in  Figs.  12  and  13  for  k  =  2. 
In  each  case,  we  present  the  results  using  two  methods  of  refinement:  uniform 
refinement  on  one  hand,  refinement  based  on  the  energy  norm  on  the  other  hand. 


c  =  0,  k  =  1 ,  unit.  ref.  c  =  0,  k  =  1 ,  unit.  ref. 


Figure  10:  Global  effectivity  index  for  the  upper  (left)  and  lower  (right)  bounds  for 
the  case  k  =  1  and  c  =  0  (top),  c  =  1  (bottom)  in  the  case  of  uniform  refinement. 


We  remark  that  the  results  in  2D  are  not  the  same  as  in  ID.  First,  the  effectivity 
indices  are  quite  different  for  p  =  1,  2  or  3,  contrary  to  what  we  obtained  in  ID. 
However,  the  effectivity  indices  for  the  upper  bound  tend  to  the  same  values  as 
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c  =  0,  k  =  1 ,  energy  norm  based  ref. 


c  =  0,  k  =  1 ,  energy  norm  based  ref. 


ndof 


Figure  11:  Global  effectivity  index  for  the  upper  (left)  and  lower  (right)  bounds  for 
the  case  k  =  1  and  c  =  0  (top),  c  =  1  (bottom)  in  the  case  of  refinement  based  on 
the  energy  norm. 


in  ID  when  p  is  increased,  that  is  1.22  for  k  =  1  and  1.07  for  k  =  2.  This  is 
an  interesting  point,  because  it  suggests  that  these  values  may  not  depend  on  the 
geometrical  dimension,  and  could  be  the  same  in  3D.  Another  interesting  remark  is 
the  performance  of  the  method  regarding  the  lower  and  upper  bounds  when  we  use 
a  refinement  based  on  the  energy  norm.  In  that  case,  the  local  problems  saturate 
very  quickly  with  respect  to  p. 


8.2  Comparison  with  other  error  estimators 

In  this  section,  we  compare  the  SRM  estimator  with  the  Element  Residual  Method 
with  (EqRM)  or  without  (ERM)  equilibrated  fluxes  and  the  Zienkiewicz  and  Zhu 
estimator  (ZZ)  based  on  a  patch-recovery  method  of  the  gradient  of  Uh.  Note  that 
the  ZZ  estimator  provides  for  an  error  estimate  only,  which  in  the  following,  is 
compared  with  the  upper  bound  estimates  obtained  from  the  SRM,  ERM,  and 
EqRM  estimators.  We  use  again  Problem  (45)  with  c  =  0.  We  choose  k  =  1  for 
the  degree  of  the  finite  element  solution  and  take  p  =  1  for  the  extra  degree  used  to 
compute  the  SRM,  ERM,  and  EqRM  estimators.  The  results  are  shown  in  Fig.  14 
on  a  sequence  of  uniformly  refined  meshes. 
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c  =  0,  k  =  2,  unit.  ref. 


ndof 


c  =  0,  k  =  2,  unit.  ref. 


ndof 


Figure  12:  Global  effectivity  index  for  the  upper  (left)  and  lower  (right)  bounds  for 
the  case  k  =  2  and  c  =  0  (top),  c  =  1  (bottom)  in  the  case  of  uniform  refinement. 


We  observe  that  the  SRM  estimator  actually  provides  for  an  upper  bound  on  the 
error  even  with  a  small  number  of  degrees  of  freedom.  Moreover  the  ZZ  estimator  is 
not  as  accurate  as  the  other  methods,  and  emphasizes  the  idea  that  this  estimator 
is  not  reliable  for  particular  problems.  Regarding  the  lower  bounds,  the  results  are 
comparable,  although  SRM  seems  to  provide  better  effectivity  indices  with  only  a 
few  degrees  of  freedom. 

8.3  Examples  of  mesh  adaptation 

The  four  methods  we  have  compared  in  the  previous  section  gave  different  estimators 
of  the  error  in  the  energy  norm.  We  now  study  the  element-wise  distribution  of  the 
error  given  by  SRM  (upper  bound),  SRM  (lower  bound)  and  ERM,  by  using  these 
error  estimators  in  the  adaptation  process  as  described  in  Section  6.  We  again 
consider  Problem  (45)  with  c  =  0,  k  =  1  for  the  degree  of  Uh,  and  p  =  1  for  the 
extra  degree  used  in  SRM  and  ERM.  Furthermore,  we  choose  Cadap  =  0.5  to  drive 
the  adaptive  process. 

The  adapted  meshes  obtained  using  the  local  contributions  to  the  global  error  are 
shown  in  Fig.  15.  We  start  with  the  same  initial  mesh  of  9  uniform  squared  elements, 
and  give  the  results  after  2,  4,  6  and  8  iterations  (from  top  to  bottom). 
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c  =  0,  k  =  2,  energy  norm  based  ref. 


c  =  0,  k  =  2,  energy  norm  based  ref. 
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Figure  13:  Global  effectivity  index  for  the  upper  (left)  and  lower  (right)  bounds  for 
the  case  k  =  2  and  c  =  0  (top),  c  =  1  (bottom)  in  the  case  of  refinement  based  on 
the  energy  norm. 


c  =  0,  k  =  1 ,  unit.  ref. 


c  =  0,  k  =  1 ,  unif.  ref. 
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Figure  14:  Global  effectivity  index  for  the  upper  (left)  and  lower  (right)  bounds  for 
the  SRM,  ERM,  EqRM,  and  ZZ  error  estimators. 


We  observe  that  the  number  of  refined  elements,  after  the  same  number  of  itera¬ 
tions,  is  larger  when  we  use  SRM  (upper  bound)  than  when  we  use  the  two  other 
methods.  This  is  due  to  the  fact  that  the  local  contributions  using  SRM  for  the 
upper  bound  are  slightly  larger  than  those  obtained  by  SRM  for  the  lower  bound 
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Figure  15:  Adapted  meshes  based  on  the  SRM  upper  bound  estimator  (left),  SRM 
lower  bound  estimator  (middle),  and  ERM  estimator  (right). 

or  ERM.  In  Fig.  16,  we  compare  the  convergence  rate  of  the  adaptation  algorithm 
based  on  the  error  estimates  provided  either  by  SRM  or  ERM.  In  both  cases,  the 
upper  bound  estimate  has  been  considered.  We  observe  that  for  a  given  number  of 
degrees  of  freedom,  the  exact  error  introduced  by  the  three  adaptation  algorithms 
is  comparable.  Figure  16  also  shows  the  different  behavior  of  the  methods  with 
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respect  to  the  number  of  elements  refined  at  each  iteration.  For  completeness,  we 
have  added  the  results  for  SRM  with  p  =  0. 


c  =  0,  k  =  1 ,  energy  norm  based  ref. 


ndof 


c  =  0,  k  =  1 ,  energy  norm  based  ref. 


Figure  16:  Graphs  showing  the  exact  error  versus  the  number  of  degrees  of  freedom 
(left)  and  number  of  degrees  of  freedom  versus  the  number  of  iterations  (right),  for 
SRM  (p  =  0  and  p  =  1)  and  ERM. 


9  Conclusions  and  future  work 

In  this  work,  we  have  investigated  a  subdomain-based  residual  method  for  a  pos¬ 
teriori  error  estimation  inspired  from  the  work  of  Carstensen  and  Funken  [6]  and 
Morin,  Nochetto,  and  Siebert  [12].  However,  our  work  is  different  in  several  aspects: 
first,  we  show  that  it  is  not  necessary  to  compute  flux  jumps  at  the  element  in¬ 
terfaces  for  the  calculation  of  the  residual.  Secondly,  we  have  tested  the  estimator 
on  quadrilateral  elements  rather  than  triangles  and  provided  a  complete  proof  of 
the  weighted  Poincare  inequality  for  this  case  (see  the  Appendix).  Moreover,  we 
have  proposed  a  method  to  recover  a  lower  bound  on  the  error  by  a  simple  post¬ 
processing  of  the  error  estimator.  We  also  show  that  reasonable  error  estimates  can 
be  obtained  by  using  polynomial  test  functions  of  same  degree  as  the  ones  used  to 
compute  the  finite  element  solution.  This  is  actually  an  important  consideration  as 
the  technology  can  then  be  implemented  in  any  existing  finite  element  codes  with¬ 
out  having  to  introduce  new  polynomial  shape  functions  of  higher  degrees,  in  other 
words,  without  having  to  drastically  modify  the  existing  data  structure.  Finally, 
the  efficiency  of  the  upper  bound  and  lower  bound  estimators  is  demonstrated  on 
ID  and  2D  elliptic  problems  and  compared  to  other  existing  methods  {i.e  the  ZZ 
method,  the  element  residual  method  and  the  equilibrated  residual  method).  The 
subdomain-based  residual  method  provides  reliable  and  robust  error  estimates  at 
reasonable  cost  for  these  test  problems.  However,  additional  numerical  experiments 
should  be  performed  to  draw  further  conclusions  about  the  quality  of  this  error 
estimator  for  broader  classes  of  problems.  For  example,  it  would  be  interesting  to 
verify  how  the  error  estimator  behaves  for  large  values  of  c  or  more  generally  in  the 
presence  of  boundary  layers. 
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Figure  17:  On  the  left,  boundary  patch;  on  the  right,  reference  element  Q  =  [0, 1]  x 
[0, 1]  divided  into  two  triangles  U  and  L. 


Nevertheless,  this  new  approach  presents  certain  attractive  features  that  need  to 
be  investigated  in  more  detail.  In  particular,  the  method  seems  suitable  for  3D 
applications  as  it  appears  to  be  much  easier  to  implement  than  the  equilibrated 
residual  method,  especially  if  the  mesh  contains  hanging  nodes.  However,  a  question 
which  needs  to  be  addressed  is  whether  the  weighted  Poincare  inequality  holds 
for  hexahedral  or  brick  elements.  Another  issue  is  whether  the  method  can  be 
easily  applied  to  linear  elasticity  and  Stokes  problems.  Additional  investigations 
on  mesh  adaptation  (refinement /derefinement)  are  advisable  as  we  have  seen  that 
for  this  method  element-wise  or  patch-wise  contributions  to  the  error  estimates  can 
be  computed.  Finally,  it  would  certainly  be  interesting  to  check  how  the  global 
error  estimator  can  be  used  for  a  posteriori  estimation  of  the  error  in  ”  quantities  of 
interest” . 


Appendix 


In  this  Appendix  we  present  the  proof  of  Lemma  2  in  the  case  of  a  2D  mesh  of 
quadrilaterals.  Two  different  proofs  are  provided,  one  for  interior  patches  and  the 
other  for  boundary  patches. 

Boundary  patch 

Let  us  consider  a  patch  uj1  on  the  boundary  (see  Fig.  17),  with  weighting  function 
4>i  that  equals  one  in  xl  and  vanishes  on  diOi  \  <9P,  and  a  function  /  vanishing  on 
<9P.  Each  of  the  two  quadrilaterals  Q±  and  Q 2  can  be  mapped  onto  the  reference 
square  Q  =  [0, 1]  x  [0, 1]  by  a  bilinear  transformation,  the  node  aq  being  mapped  to 
(1,1).  The  weight  (j>i  is  then  transformed  into  </>  =  xy.  Then,  by  a  scaling  and  a 
density  argument,  it  is  sufficient  to  prove  the  inequality 

[  f2(x,y)dxdy  <C  [  \V  f(x,y)\2xy  dxdy  (47) 

JQ  JQ 


for  any  smooth  function  /  vanishing  either  on  y  =  1  or  x  =  1. 
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Let  us  consider  the  case  where  /  =  0  at  y  =  1  (the  other  case  can  be  proved  in  a 
similar  manner).  We  establish  the  proof  of  (47)  in  two  steps.  The  basic  idea  of  this 
proof,  which  is  to  integrate  by  parts  along  a  line  y  —  x  =  constant,  is  borrowed  from 
the  paper  by  Carstensen  et  al.  [6,  Lemma  5.3]. 


Step  1.  We  first  prove  the  inequality 


>u 


f2(x,y)dxdy<C1  |  V/(x,  y)\2xy  dxdy 


>u 


where  U  is  the  upper  triangle  shown  in  Fig.  17. 


(48) 


By  introducing  the  change  of  variables  x'  =  x  and  y'  =  y  —  x ,  which  transforms  the 
triangle  U  =  {(0, 0),  (1, 1),  (0, 1)}  into  the  triangle  U'  =  {(0,  0),  (1,0),  (0, 1)},  and  by 
integrating  by  parts  along  each  line  y'  =  constant,  we  can  write 


lu 


f2(x,y)dxdy=  [  dy'f  f2{x',y')dx' 

Jo  Jo 

fi 


=  /  dy 


f  {x  ,y  )x‘ 


l~y'  fX~V\  rdf  ; 

-  /  2finx  dx 

o  Jo  dx 


=  -  f  2f(BJ-  +  ^\Xd,iy. 

Ju  \ox  ay 


The  boundary  term  vanishes  at  both  ends.  In  particular,  the  line  x'  =  1  —  y' 
corresponds  to  y  =  1  and  f\y=1  =  0.  We  finally  have 

f2  dxdy  <  ^  j^f2  dxdy  +  2  x2  dxdy 

<  \  f  f  dxdy  +  4  [  \X7  f\2xy  dxdy 

z  Ju  Ju 

where,  in  the  last  inequality,  we  have  used  the  discrete  Cauchy  inequality  lai|  5s 

y/n  af)1^2  and  the  fact  that  x2  <  xy,  V  (x,  y)  €  U.  This  completes  the  proof 

of  inequality  (48)  with  a  constant  C\  =  8. 


Step  2.  We  now  prove  the  inequality 


/  f2(x,y)dxdy  <  C2  |V  f{x,y)\2xy  dxdy 
Jl  Jq 

where  L  is  the  lower  triangle  shown  in  Fig.  17.  Indeed  we  have 


rl  rx 

f2(x,  y)  dxdy  =  dx  f2(x,y)dy 

Jo  Jo 

f2{x,  y)y 


=  dx 


-  [X2fdMdy 
Jo  oy 


(49) 


<  [  f2(x,x)xdx  +  i  f  f2  dxdy  +  2  [  y2 dxdy , 

Jo  2  Jl  JlKovJ 


28 


Prudhomme,  Nobile,  Chamoin  and  Oden 


thus  leading  to  the  inequality 


f2(x,y)dxdy<  2  f  f2(x,  x)x  dx  +  4  f  |V  f(x,y)\2xy  dxdy  (50) 

Jo  Jl 


where,  again,  we  have  exploited  the  fact  that  y2  <  xy,  V  ( x ,  y)  £  L. 

To  estimate  the  boundary  term  f2(x,  x)x  dx  appearing  in  the  previous  inequality, 
we  proceed  as  follows: 


f2(x1y)dxdy 


f2(x,y)dx 


f2{x, y)x 


of9f  ,  ' 

2t-—xdx 
J  dx 


f2{y,y)ydy 


f  2 f^-xdxdy. 

Ju  ox 


Hence, 


^  f2{y,  y)y  dy  =  f2{x,  y)  dxdy  +  £  2f^x  dxdy 


iu 


<2  f  f2(x,  y)  dxdy  +  [  \Vf\2xydxdy 
Ju  Ju 

<(2Ci  +  l)  f  \V f\2xy dxdy 
Ju 


(51) 


where,  in  the  last  inequality,  we  have  employed  relation  (48).  Finally,  by  inserting 
(51)  in  (50),  we  obtain  (49)  with  C*2  =  2(2C*i  +  1),  and,  by  summing  (48)  and  (49), 
we  finally  establish  the  statement  claimed  earlier. 


Interior  patch 

Let  us  now  consider  an  interior  patch  u>i  of  quadrilaterals  (see  Fig.  18).  In  this  case, 
the  weight  function  (f>i  equals  one  at  Xi  and  vanishes  on  duii,  while  the  function  /  is 
such  that  f4>i  =  0.  The  patch  Ui  can  be  mapped  onto  the  reference  patch  Qj  = 
[— 1, 1] 2 ,  shown  in  Fig.  18,  by  a  piecewise  bilinear  transformation  \  u  —>  LOi.  The 
node  Xi  is  mapped  to  (0,  0)  and  the  weight  <f>i  is  transformed  as  =  (1  —  |x|)(l  —  |y|). 
Then,  by  a  scaling  and  a  density  argument,  it  is  sufficient  to  prove  the  inequality 


/  -  -  /  f4>Ji  duj 
y  Jco 


L2(u) 


<  C  /  |V/|2<^do;,  d=  <JJiduj, 


(52) 


for  any  smooth  function  /,  with  a  constant  C  that  may  depend  on  the  regularity  Ki 
of  the  patch  c Oi  but  is  otherwise  independent  of  hi.  We  have  denoted  Jj  =  det(V^i). 
Actually,  the  function  g  =  f  —  ^  f(f>Ji  satisfies  the  condition  gfa  =  J'^  g<l>Jt  =  0 
and  (52)  reads  also 
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Figure  18:  Interior  patch  of  quadrilaterals  (left)  and  patch  of  reference  (right).  Each 
square  in  the  patch  of  reference  has  been  divided  in  two  triangles. 


We  will  prove  the  inequalities 


/  -  -  [  f&Ji  duJ 


L2{Tj) 


<Cj 


f  |V/|‘ 

J  UJ 


<l>du,  V  j  =  1, ...  ,8,  (53) 


from  which  (52)  is  easily  obtained  by  summing  over  the  indices  j.  The  proof  is 
carried  out  only  for  the  triangle  Tf,  the  other  cases  being  obtained  in  a  similar 
manner. 

Let  us  introduce  the  change  of  variables  x'  =  1  —  x  and  y'  =  x  —  y,  which  trans¬ 
forms  the  triangle  T\  =  {(0, 0),  (1,0),  (1, 1)}  into  {(1,0),  (0, 1),  (0, 0)}.  Then,  for  any 
smooth  function  g  on  T\,  we  have 


IT i 


/•I  r-L-y 

g2(x,y)dxdy  =  dy'l  g2(x'  ,y)  dx'dy' 

Jo  Jo 

=  f  dy' 

Jo 


i—y'  rl~y'  da 
o  2aMX'dx 


< 


g\x',y')x ' 

g2( 1  -  y',  y')(l  -  y')  dy'  +  2  g  (1  -  x)  dxdy 

[  g2(x ,  0)(1  —  x)  dx 

Jo 

+  \  j  g2  dxdy +  4  [  | V^|2(l  -  x)(l  -  y)  dxdy. 

1  Jt-l  Jt i 


By  taking  g  =  f  —  ^  dco,  we  obtain  the  inequality 


1 


/ - /  dto 

d  J  co 


L2(Ti) 


fS’-lL 


<  2  /  / - /  f4>Ji  dw  |  (1  —  x)  dx  +  8  /  |V/|  i j)duj 


li/=0 


’T\ 
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We  now  need  to  estimate  the  boundary  term  appearing  on  the  right-hand  side  of 
the  previous  inequality. 


f  -  -  [  f4>Ji  dcu)  |  (1  -  x)  dx 

/0  V  d  Jul  J  I  y=0 


1 


/  0  I1 


(f(x,  0)  —  f(t,  s ))  (j)(t,  s)Ji(t,  s)  dtds  |  (1  —  x)  dx 
/  (1  —  x)  dx 

Jo 

(/  -^-(x,v)dv  +  J  s)  du^j  s)Ji(t,  s)  dtds 


<  —  max|Jj|2(/i  +  /2)  <  C(Ii  +  I2) 

H z  ^ 


where 


Ii 

/2 


f\l-x) 

Jo 

f\l-x) 

Jo 


df 

—(x,v)<t)(t,s)dv 


2 


dx, 


df 

—  (u,s)<j>(t,s)du 


2 


dx. 


Observe  that  the  constant  C  depends  only  on  the  regularity  n;  of  the  patch  u;,  since 
\Ji\<  Clif  and  n  >  Cpj  f. 

We  separately  analyze  the  two  terms  I±  and  I2.  We  start  with  I2\ 


I2  =  /  (1  —  x)dx 


i>i: 


d f 

du 


( u ,  s )  y /  <f(u,  s)  yi  1  =  du 
f{u,s) 


ds 


<  /  (1  —  x)dx 
Jo 


L 


df ,  x 


s)  duds 


r*l  /-X  i2 


(t,s) 

I- 1  Jt  (j>(u,  s) 


duds 


dt 
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that  is, 


(1-IOT-I.DJ” 

1  ~U 


In  2 


<c\  I  \\7 f\24>du 

J  LJ 

where  the  constant  C\  is  given  by 


C\  =  f  (1  —  x)  dx 

Jo  [J-i 

r  i  r  r  i 

<  (1  —  x)  dx 

Jo 


/  r1 

r  I,* 

i  -i 

\  2 

(  /  (1  -  \s\ )ds 

\J- 1 

1 

— i 

J 

(1  -  |i|)  [  log(l  -  |x|)  log(l  -  |t|)]2  dt 


<  +oo. 


We  now  consider  the  term  I\\ 


I\  =  ^  dx  j  j  dt  J  ds 


y^(x,v)\J '  dv 


x,v) 


<  [  dx 

Jo 


/>/>/, 


Bf 

dv 


(x,v) 


(j)(x,  v )  dv 


rO  12 


<r(M)(  1  -  |x|) 

<i>{x,v) 


dv 


i~\  2 

2 


r*l 


r*l 


-L  * \L 


01 

dv 


(x,v) 


j  dt  j  ds 


<j){x,  v)  dv  J 

0  (i  -  l*l)2(i  -  M)2 
(i  -  H) 


dv 


In  2 
2 


<C2  \V  f\2(j)duj 

J  LJ 

where  the  constant  C2  is  given  by 


C2  = 


j  (1  —  \t\)  dt  J  (1  —  |s|)ds 


(1  -  |a|)  [-  log(l  -  |a|)] 2  ds  <  +oo 


1  -  \v\ 


■  dv 


In  2 

2 
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Observe  that  the  constant  C2  is  bounded  thanks  to  the  presence  of  the  weight  (1 — x). 
We  finally  obtain 

f1  (f-  -  [  fHuj)  I  (1  -  x)  dx  <  C(C 1  +  C2)  [  \V f\2cj)dco. 

JO  V  M  Ju)  J  I y=0  J Cb 

Then,  inequality  (53)  is  satisfied  on  the  triangle  Tj.  An  analogous  proof  holds  for 
the  other  triangles,  thus  leading  to  the  final  result. 

Acknowledgement.  The  support  of  this  work  through  ONR  grant  N00014-95-0401 
is  gratefully  acknowledged.  The  authors  wish  also  to  acknowledge  Prof.  Babuska 
for  the  fruitful  discussion  on  the  subject. 


References 

[1]  M.  Ainsworth  and  J.T.  Oden,  A  Posteriori  Error  Estimation  in  Finite  Element 
Analysis,  John  Wiley  and  Sons,  2000. 

[2]  M.  Ainsworth  and  J.T.  Oden,  A  unified  approach  to  a  a  posteriori  error  esti¬ 
mation  based  on  element  residual  methods,  Numer.  Math.,  65,  23  (1993). 

[3]  I.  Babuska  and  W.C.  Rheinboldt,  A  posteriori  error  estimates  for  the  finite 
element  method,  Internat.  J.  Numer.  Methods  Engrg.,  12,  1597  (1978). 

[4]  I.  Babuska,  T.  Strouboulis  and  S.K.  Gangaraj,  Guaranteed  computable 
bounds  for  the  exact  error  in  the  finite  element  solution.  Part  I:  One- 
dimensional  model  problem,  Comput.  Methods  Appl.  Mech.  Engrg.,  176,  51 
(1999). 

[5]  R.E.  Bank  and  A.  Weiser,  Some  a  posteriori  error  estimators  for  elliptic  partial 
differential  equations,  Math.  Comp.,  44,  283  (1985). 

[6]  C.  Carstensen  and  S.A.  Funken,  Fully  reliable  localized  error  control  in  the 
FEM,  SIAM  Journal  of  Scientific  Computing,  21,  1465  (2000). 

[7]  P.G.  Ciarlet,  The  Finite  Element  Method  for  Elliptic  Problems,  North- 
Holland,  1978. 

[8]  D.K.  Datta,  Computer  Analysis  of  Error  Estimation  in  Finite  Element  Com¬ 
putations  for  Elliptic  and  Parabolic  Problems,  Ph.D.  Dissertation,  Texas  A&M 
University  (2001). 

[9]  L.  Demkowicz,  J.T.  Oden,  and  T.  Strouboulis,  Adaptative  finite  elements  for 
flow  problems  with  moving  boundaries.  Part  1:  Variational  principles  and 
a  posteriori  error  estimates,  Comput.  Methods  Appl.  Mech.  Engrg.,  46,  217 
(1984). 


Subdomain-Based  Residual  Method 


33 


[10]  A.  Ern  and  J.-L.  Guermond,  Elements  Finis:  Theorie,  Applications,  Mise  en 
Oeuvre,  Springer- Verlag,  Berlin,  2002. 

[11]  P.  Ladeveze  and  D.  Leguillon,  Error  estimate  procedure  in  the  finite  element 
method  and  applications,  SIAM  J.  Numer.  Anal,  20,  485  (1983). 

[12]  P.  Morin,  R.H.  Nochetto,  and  K.G.  Siebert,  Local  problems  on  stars:  a  pos¬ 
teriori  error  estimators,  convergence,  and  performance,  Math.  Comp.,  (to  ap¬ 
pear). 

[13]  T.  Strouboulis,  I.  Babuska,  and  S.K.  Gangaraj,  Guaranteed  computable 
bounds  for  the  exact  error  in  the  finite  element  solution.  Part  II:  bounds  for 
the  energy  norm  of  the  error  in  two  dimensions,  Internat.  J.  Numer.  Methods 
Engrg.,  47,  427  (2000). 

[14]  O.C.  Zienkiewicz  and  J.Z.  Zhu,  The  superconvergent  patch  recovery  and  a 
posteriori  e rror  estimates.  Part  1:  The  recovery  technique,  Internat.  J.  Numer. 
Methods  Engrg.,  33,  1331  (1992). 

[15]  O.C.  Zienkiewicz  and  J.Z.  Zhu,  The  superconvergent  patch  recovery  and  a 
posteriori  error  estimates.  Part  2:  Error  estimates  and  adaptivity,  Internat.  J. 
Numer.  Methods  Engrg.,  33,  1365  (1992). 


